Analysis of vibration impact on stability of dewetting thin liquid film 
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Dynamics of a thin dewetting liquid film on a vertically oscillating substrate is considered. We 
assume moderate vibration frequency and large (compared to the mean film thickness) vibration 
amplitude. Using the lubrication approximation and the averaging method, we formulate the coupled 
sets of equations governing the pulsatile and the averaged fluid flows in the film, and then derive the 
nonlinear amplitude equation for the averaged film thickness. We show that there exists a window 
in the frequency-amplitude domain where the parametric and shear-flow instabilities of the pulsatile 
flow do not emerge. As a consequence, in this window the averaged description is reasonable and 
the amplitude equation holds. 

The linear and nonlinear analyses of the amplitude equation and the numerical computations 
show that such vibration stabilizes the film against dewetting and rupture. 

PACS numbers: 47.15.gm, 47.20.Ma, 68.08.Bc 



I. INTRODUCTION 

Studies of stability, dewetting and rupture of thin 
liquid films on solid substrates are of great impor- 
tance for micro- and nanotcchnologies. Since the first 
works appeared in 1960s and 1970s, these subjects con- 
tinue to attract enormous attention, see for instance 
Refs. [1 0, i 1, i, S 0, |, i, [p. Refs. and 
are the reviews of recent results in this extremely diverse 
field. 

The subject of this paper is the general theoretical in- 
vestigation of the impact of the vibration on the sta- 
bility of a thin dewetting liquid film. Our interest in 
studying vibration impacts stems from the large body 
of works in fluid mechanics of macroscopic fluid layers, 

EE EH (the experiment) and 
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20) [2l| (the theory and numerical mod- 



including Refs. 
Refs. [HEIM 
eling), where the vibration is shown to drastically affect 
the stability characteristics and the dynamics of fluid sur- 
faces and interfaces. 

In this paper we assume that the external influences on 
the film are the vibration, the gravity, and the long-range 
molecular attraction by the planar substrate, typified by 
van der Waals forces. Other effects (such as, for instance, 
the thermocapillarity and the evaporation) can be easily 
included by simply adding the corresponding terms to 
the final evolution equation for the film thickness. 

The theory we develop is based on the standard long- 
wave lubrication approximation, as discussed by Oron, 
Davis, and Bankoff [111 ], and on the time-averaging 
method. The general discussion of the averaging meth- 
ods can be found in Refs. (22l , l23l . Izij . The key idea is the 
separation of the dynamics onto fast pulsations and slow 
relaxation processes. This approach works well when the 
vibration frequency is high in a certain sense, i.e. when 
there exists a large difference in the characteristic times 
(such as the viscous relaxation time and the vibration 
period, see below). 



The first transparent explanation of such separation of 
the time scales was given by Kapitza in his pioneering 
work on a pendulum with an oscillating point of sup- 
port . The paper by Blekhman [26| contains other 
examples in mechanics. Many examples of the successful 
application of the averaging method can be found in ther- 
mal vibrational convection (l7j . dynamics of inclusions 
in fluids [H, H3, HH , dynamics of granular materials [2!| , 
motion of disperse fluids [3(| , and filtration of inclusions 
in porous media [3l| . 

Interestingly, the vibration of the solid plate (on which 
the fluid system is located) often is capable of complete 
suppression of instabilities. For example, Wolf [1J] ex- 
perimentally investigated the damping of the Rayleigh- 
Taylor instability in the horizontal two-layer system by a 
vertical high frequency vibration. The theoretical anal- 
ysis of this situation (in the linear approximation) by 
the averaging technique was performed by Cherepanov 
[32| (summary of this paper can be found in Ref. [Hj]). 
The longwave instability in this system was also ana- 
lyzed by Lapuerta, Mancebo, and Vega [l§|. After the 
analysis of the linear longwave instability at the mod- 
erate vibration frequency, they proceed to the averaged 
description at high frequency. The generalization of the 
latter analysis to the nonisothermal situation was devel- 
oped by Thiele, Vega, and Knobloch [20]. They account 
for the Marangoni effect and perform a detailed investi- 
gation of the corresponding amplitude equation. 

Another widespread vibration-induced phenomenon is 
the parametric excitation, which emerges when the fre- 
quency of the vibration is comparable to one of the eigen- 
frequencies of the system (for instance, to the frequency 
of the capillary-gravity waves). Faraday was first to ob- 
serve parametric waves on the surface of vertically oscil- 
lating horizontal layer (33j . Linear and nonlinear analyses 
of parametric instability were performed, for instance, by 
Benjamin & Ursell [34L Kumar & Tuckerman [3^ . Lyu- 
bimov & Cherepanov [36| . and by Mancebo and Vega 
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[37l |. To the best of our knowledge the latter paper is 
the most detailed study to-date of the linear aspects of 
Faraday instability. 

It must be noted that the situations termed "the aver- 
aged motion" and "the parametric instability" are often 
closely connected, although they operate within different 
intervals of the vibration frequencies. Indeed, in the stud- 
ies of the averaged dynamics one has to ensure stability 
of the pulsatile motion (periodic in time) . Most fluid sys- 
tems have an eigenfrequency spectrum unbounded from 
above and thus the eigenfrequency is an increasing func- 
tion of the mode number. Thus, even the high-frequency 
vibration is capable of parametric excitation of the higher 
modes. For the pulsatile motion to be stable, a window 
of parameters such as the amplitude and the frequency 
of the vibration must be chosen, where the parametric 
instability does not emerge. 

It is also worth noting that most papers [13, EH H3| , 
where the averaging method is employed, deal with the 
vibration of "inviscid" frequency, i.e. the vibration pe- 
riod T p — 2it/u; (here to is the dimensional frequency) is 
assumed small compared to the characteristic time of vis- 
cous relaxation, t v — H 2 /v (here v is the kinematic vis- 
cosity and Hq is the mean fluid layer thickness) . It is clear 
that this assumption is quite reasonable for macroscopic 
layers, but for thin films (of thickness 100— 1000 A) it re- 
quires extremely large frequencies, 100 MHz and higher. 

Now, we make a very important point, as follows. A 
thin film allows for averaged description even when the 
viscosity is large [i.e., t v = O (T p ) and even t v <C T p j: 
one needs only to assume that the period of the vibra- 
tion is small compared to the characteristic time of the 
film evolution, T\ = O (fc 2 ) (here k is the wavenumber). 
Due to the lubrication approximation, t; S> t v . Thus 
the condition t; 3> T p is much milder than the usual in- 
viscid approximation t v 3> T p . Therefore, the averaging 
procedure can be applied even to ultrathin films. 

To the best of our knowledge the only paper devel- 
oping similar approximation is Ref. [19j], where the lin- 
ear stability problem is studied for "moderate" ( "finite" ) 
non-dimensional frequencies lot v = 0(1). However, by 
assuming the amplitude of the vibration "finite" (which 
means that it is of the order of the fluid layer thickness), 
the authors obtain that the impact of the vibration at 
moderate frequency is small. Thus, they focus on the 
high-frequency, small- amplitude case lut v 3> 1. 

In this paper we assume large vibration amplitude and 
develop the nonlinear amplitude equation for the thick- 
ness of the film. 

The outline of the paper is as follows. The mathemat- 
ical formulation of the problem is presented in Sec. HT1 
We state and discuss assumptions and governing equa- 
tions. In particular, the separation of the processes onto 
fast pulsatile and slow averaged motions is performed in 
Sec. Ill Bl The pulsatile problem is solved in Sec. IIIII The 
amplitude equation [Eq. (|36[) ] governing the averaged dy- 
namics of the film thickness is obtained in Sec. IIVI This 
equation is the main result of this paper. It can be used to 
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FIG. 1: Problem geometry. 



study impacts of the vertical vibration (in the frequency 
range for which the averaged description is applicable) 
on the dynamics of a film, in the presence of the sur- 
face tension and wetting interactions with the substrate. 
Two limiting cases of "low" and "high" vibration fre- 
quency are analyzed in Sec. [V] These cases correspond 
to the different ratios of the inertial and viscous forces 
in the oscillatory motion. (The viscosity dominates at 
low frequencies, while the inertia force dominates at high 
frequencies.) In Sec. I VII the 3D generalization of the 
theory is presented. Conditions of parametric instability 
of the oscillatory motion are analyzed in Sec. IVIH where 
the Faraday instability and the shear flow instability are 
discussed. We show that for any admissible vibration fre- 
quency there exists a finite range of vibration amplitudes 
for which such instabilities are not present. In Sec. IVIIIl 
we address the averaged behavior of the system within 
the framework of the obtained amplitude equation. Re- 
sults of linear and weakly nonlinear analyses of the equi- 
librium state with the flat surface are presented, as well 
as the results of direct numerical simulations. In partic- 
ular, we show that the vibration influence is stabilizing, 
i.e. it can delay or completely suppress the film rupture 
by intermolecular attractive forces. Finally, Sec. lIXI sum- 
marizes the results. 



II. FORMULATION OF THE PROBLEM 
A. Governing equations 

We consider a three-dimensional (3D), laterally un- 
bounded thin liquid film of unperturbed thickness Ho 
on a planar, horizontal substrate. The Cartesian refer- 
ence frame is chosen such that the x— and y— axes are 
in the substrate plane and the z— axis is normal to the 
substrate (Fig. [TJ. 

The substrate-film system is subjected to the vertical 
harmonic vibration of the amplitude b and the frequency 
ui. Thus, in the reference frame of the substrate the ac- 
celeration of gravity is modulated, 

g(t) = go + buj 2 cos tot. (1) 

Since H is small, the intermolecular interaction of the 
film surface with the substrate has to be taken into ac- 
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count. Thus, we prescribe the potential energy (f>(H) to 
the unit length of the film layer. In this paper we consider 
only the van der Waals attractive potential, 



(H) 



A' 



(2) 



where A' > is the Hamaker constant [38|]. The model 
presented in this paper can be readily extended to incor- 
porate other models of wetting interactions - it is only 
necessary to replace 4>(H) with an appropriate function. 

We scale the time, the length, the velocity and the pres- 
sure by Hq/v, H , v/Hq, p(v/H ) 2 , respectively (here, 
v is the kinematic viscosity and p is the density of the 
liquid). Then, the liquid motion is governed by the fol- 
lowing non-dimensional problem: 

Vv = 0, (3a) 
vj+v-Vv = -Vp + V 2 v 

- (Go + Bfl 2 cos fit) e z , (3b) 

v = at z = 0, (4a) 
ht + v • V/i = w, 
(p- 4>- CaK)n = n • T at z = h(x,y,t). (4b) 

Here, v is the fluid velocity, w is its component normal 
to the substrate, p is the pressure in the liquid, T is 
the viscous stress tensor, h — H/Hq is the dimension- 
less thickness of the layer, e z is the unit vector directed 



along the z— axis, n = (e z — Vh) /y 1 + (V/i) is the 
normal unit vector to the free surface, K = V ■ n is 
the mean curvature of the free surface, cf) — A/h 3 [where 
A = A' I (6ir pv 2 Ho) is the non-dimensional Hamaker con- 
stant], Ca = <jHq/ '{pv 2 ) is the capillary number (where 
a is the surface tension), Go = gaH^/v 2 is the Galileo 
number, B = b/Ho is the non-dimensional amplitude, 
and fl = loHq jv is the non-dimensional frequency. 

We consider the nonlinear evolution of the large-scale 
perturbations. As it is usually done, we introduce small 
parameter e, which is of the order of the ratio of the mean 
thickness Hq to the perturbation wavelength, i.e. e <C 1 
for long waves. 

Below for the sake of simplicity and for more trans- 
parent presentation of ideas, we consider the 2D model, 
assuming that v = ue x + we Zl where e x is the unit vector 
directed along the x— axis, and all fields are independent 
of y. The theory extends trivially in three dimensions at 
the replacement of the x-derivative by the 2D gradient. 
We derive the 3D analogue in Sec. IVII 

Next we introduce the conventional stretched coordi- 
nates and the time: 



X = ex, Z = z, T = e 2 t, r = fit, 



(5) 



such that d/dt — e 2 d/dT + 0,8/ dr. Now we rescale the 
velocity components as follows: 



To this end, both the pressure p and the surface position 
h remain unsealed. 

Throughout this paper we assume that the capillary 
number is large: 



Ca = Ce" 2 . 



(7) 



This is quite realistic and widely used assumption [111 ]. 

Substituting the expansions ([5])- (J7J) into the problem 
([3]) and ([4]) we arrive at the following set of equations 
and boundary conditions: 

W z = -U x , (8a) 
nu T = -PX + UZZ 

+e 2 (U xx ~U T - UU x - WU Z ) , (8b) 
pz = -G -BQ?casT + e 2 {Wzz-O.W T ), (8c) 



U = 0, W = at Z = 0, 

U z = -e 2 (W x ~ ihxUx) , 

flh T = e 2 (W - h T -Uh x ) , 

p = <j>- Chxx 



(9a) 



-e 2 ( ~Ch 2 x hxx + 2W Z ) at Z = h. (9b) 



u = ell, w = e 2 W. 



(6) 



All terms of order e 4 have been omitted. 



B. Separation of the pulsating and averaged 
dynamics 

In this paper we consider the case of moderate vibra- 
tion frequency: fl 3> e 2 , i.e. the vibration period is small 
compared to the characteristic time of the surface evo- 
lution. This assumption makes possible the averaging 
of the dynamics of the film over the vibration period 
HH HI, 0] . The main purpose is the rigorous asymptotic 
analysis (in powers of e) which results in the set of equa- 
tions and boundary conditions, where the dependence on 
r is averaged out. This set is then used to derive the am- 
plitude equation for the thickness of the film. As for fi, 
this quantity is not an asymptotic parameter. This means 
that, generally speaking, we assume £1 neither large nor 
small. Therefore Q will enter the amplitude equation as 
a finite parameter. 

However, within the framework of the equation set 
([5]) and © we will consider three qualitatively different 
cases: 

(i) 12 <gC 1, i.e. the vibration is of "low" frequency; this 
is the quasi-Stokes approximation with week influence of 
fluid inertia. Note that in this case, due to Eqs. |5]), 
the condition e 2 <SC il must be retained to legalize the 
averaging procedure. 

(ii) = O(l). This assumption means that the vibra- 
tion period is comparable to the time of the momentum 
relaxation across the layer, Hq / v. 

(iii) S! 3> 1, i.e. the vibration is of "high" frequency; 
viscosity is negligible except for the thin boundary layer 
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near the rigid wall. As usual, the boundary layer in the 
vicinity of the free surface does not play an important 
role (see Ref. [39 ]. for example). 

Since the van der Waals interaction is important for 
films of thickness 100— 1000 A, let us estimate the typical 
values of f2. Taking Hq = 10~ 5 cm and v = 10~ 2 cm 2 s _1 
(viscosity of water), we obtain « 10~ 4 for uj/(2tt) = 
1 kHz (typical for mechanical vibrator), and f2 = 10~ 2 -^1 
for n = 0.1 -j- 10 MHz (typical for ultrasound). Therefore 
the case (i) can be easily achieved either mechanically or 
by means of ultrasound irradiation of the substrate from 
below, the case (ii) can be reached only using the ultra- 
sound, and the case (iii) seems unrealistic. Nevertheless 
we shall consider this case below, since our results, upon 
neglecting 0, can as well be applied to the description of 
macroscopic fluid layers. Besides, this limiting case was 
studied in detail by many authors and therefore it allows 
for the verification of our results. 

We represent each field / = {U,W,p,h} as the sum 
of the average part f(T) and the pulsation /(r), where 
r and T can now be termed the "fast" and the "slow" 
times, respectively [24j. Thus, we write 

U = U + U, W = W + W, (10a) 
p = p + p 1 h = h + h. (10b) 

We now assume that 

B = e- 1 b, (11) 

i.e. the amplitude of the vibration is large compared to 
the mean film thickness. This assumption seems sur- 
prising as it is customary to impose the small- amplitude 
high frequency vibration. However, large amplitudes are 
permitted when the large-scale dynamics is considered. 
Indeed, it is shown below that due to the longwave ap- 
proximation the impact of the vibration becomes non- 
negligible only at large amplitudes. Also, it will be made 
clear momentarily that in some sense the pulsatile motion 
is still small-amplitude. 

The assumption of large vibration amplitude means 
that the oscillating part of the pressure field p (which 
is forced by the inertia force -Bf2 2 cosr) is of order e _1 , 
which in turn leads to e^ 1 scaling for the velocities of the 
pulsation U and W . The pulsation of the surface height 
h obviously has the same order as w — e 2 W [see Eq. JBJ], 
i.e. it is of the order e. Therefore, it is convenient to re- 
define the pulsation parts of all fields, rewriting Eqs. (flU]) 
as follows 

U = U + e^U, W = W + e~ 1 W, (12a) 
p = p + e~ 1 p, h = h + eh, (12b) 

where U,W,p and h are O(l) quantities. 

Accounting for the initial scaling © one can conclude 
that the full components of the velocity field are: 

u = eU + U, w = e 2 W + eW, (13) 
while the pressure field and the surface deviation are 
given by Eqs. (| 12b|) . Note that the pulsations of the 



fluid velocity and the pressure are larger than their mean 
parts, while the opposite is true for the pulsation of the 
film height. Moreover, the scaling (fT3"|) means that the 
pulsation remains "small-amplitude" . Indeed, the typ- 
ical horizontal (vertical) displacement of the fluid par- 
ticle during one period is O(l) [O(e)], which is small in 
comparison with the respective characteristic lengthscale, 
0(e-i) [0(1)]. 

Substitution of the expansions (fT2|) in the equation sets 
© and © allows to separate fast pulsations from back- 
ground slow (averaged) motions. Keeping terms of zeroth 
and first orders in e, we obtain the following sets. 

(i) For the pulsations: 



W z 

VLU T 

Pz 
Uz 



-U x , 

-Px + U ZZ -e (utJx + 

—bfl 2 COS T, 

0, W = at Z = 0, 

Uzhhx 



WUz) ■ 

/ P 



-Uh x + W - 
-e(u zz h) .. 



(14a) 

(14b) 
(14c) 



x 



P : 



-e(pzh) p at Z = /<14d) 



where the subscript "p" denotes the pulsating part of the 
corresponding expression: 



f P = .f- /• 



(15) 



(It is obvious that in general the term squared with re- 
spect to pulsations contains both the averaged and the 
pulsating components, for instance cos 2 r = \ + \ cos2r.) 
(ii) For the averaged parts: 



W z 

Uzz 
U 

P 
U z 



-U x , Pz = -G , (16a) 

Px + (UU x + WUz), (16b) 

0, W = at Z = 0, (16c) 

-(Pzh) - Upzzh 2 ) + <f>(h) - Chxx, 



^Uhx - (Uh x ) + W+ (W z h) - Uu z h 2 ) x , 



-(Uzzh)- e -(U Z zzh 2 



at Z = h. 



(16d) 



In the system (fT6)) the angular brackets denote averaging 
in time r. Note that the boundary conditions at the free 
surface have been shifted at the mean position h. This 
leads to the following expansion in powers of e of the 
arbitrary field F: 

F(Z = h + eh) w F{h) + eF z (h)h + ~e 2 F zz (h)h 2 . (17) 



As we have noted, the e 2 -terms in Eqs. (fl6|) have been 
omitted. On the other hand, the third term in Eq. (fTTI) 
has to be taken into account for some fields, because it 
produces a correction of order e. For instance: 
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p(Z = ft + eft) 



e 1 [p + ep z h + -e 2 p zz h 2 ) n + <p z h 

p(Z = ft) + (p z ~h + \-tp zz h 2 

\ 1 I Z =h 



Z=h 



(18) 



r 



However, all terms cubic with respect to pulsations van- 
ish after the averaging. Moreover, in further analysis we 
will disregard the e-terms in Eqs. (fl"4|) and (fl~6|) (i.e., use 
these boundary value problems in the leading zeroth or- 
der). 

It is also important to emphasize that the zeroth order 
of the problem governing the pulsations [problem (|14|)] is 
linear. This is because all nonlinear terms (UUx, WUz, 
etc.) are small due to the longwave approximation, de- 
spite the scaling (fT3)) . 



III. PULSATILE MOTION 

Equations. (|14p become in the zeroth order in e: 

p z = -6f2 2 cosT, (19a) 

W T = ~px + U zz , W Z = ~U X , (19b) 

U = IT = at Z = 0, (19c) 
p = 0, U z = 0, 

Qh T = -Uh x + W at Z = ft. (19d) 

According to the conventional method for the solution 
of linear problems, we need to separate the general so- 
lution of the homogeneous problem from the particular 
solution of the nonhomogeneous one. The former solu- 
tion corresponds to the gravity-capillary waves damped 
by viscosity and thus it is not of interest. (In fact, the 
gravity-capillary waves are completely damped at times 
of order r. Thus the omitted solution represents fast 
relaxation of the initial conditions, which has no effect 
on slow dynamics with the characteristic time scale T.) 
The latter solution corresponds to the periodically forced 
motion, which can be represented in the complex form as 
follows: 



where a 2 — Hereafter we omit the bar over ft. 

At a — this problem coincides with the conventional 
equations for the thin film in the absence of the vibration 
[cf. Eqs. (2.22)-(2.24) in Ref. (Tj at $ = (3 = £ = n = 
C _1 = 0]. However, the term a 2 1 originating from the 
inertia of the fluid drastically complicates the solution. 
The solution of the boundary value problem (j2"Tj) is: 



q = Q(h- Z), I = ih x [l 

K = - 



i>xx 



Z 



cosa(/i — Z) 
cos ah 
sina(/i — Z) — sin aft, 



a cos ah 



h 2 x 



1 — cos aZ 



H = 



cos 2 ah 

— \hhxxf{oth) — h\ tan 2 ah\ 

- [hf(ah)h x ] x , 



where 



tan y 
V 



(22a) 



(22b) 

(22c) 
(22d) 



(23) 



Now, if we set h = 1 + £ and linearize Eqs. (|22|) with 
respect to £, then we arrive at Eqs. (2.30)-(2.32) in 
Ref. [19(. However, we have considerably expanded the 
domain of validity for this solution, as will be explained 
in Sees. |W] and Em! 

Next, we proceed to the analysis of the limiting cases 
for the pulsatile motion, i.e. Q <C 1 (low frequency) and 
f2 3> 1 (high frequency). 

In the former case the viscous term Izz dominates 
in Eq. (|21a[) and the inertial term — iQI exerts a week 
impact only. Thus the solution of the problem (f2"Tj) sim- 
plifies to: 



p = bflRc [q(X, Z) exp (ir)} , 
U = bflRe [I(X, Z) exp (ir)} , 


(20a) 
(20b) 


q = 
i = 


fi(ft- 

nh x 

2 


W = bttRe [K(X, Z) exp (ir)} , 


(20c) 




X 


1 - 


h = bRc [H(X) exp (ir)} . 


(20d) 




The set of equations and boundary conditions governing 
the amplitudes of pulsations reads: 


K = 


Vlhxx 
6 

nh 2 x 


qz = Izz +a I = q x , K z = -Ix 


(21a) 






2 


I = K = at Z = 0, 

q = 0, iH = K - Ih x , /z = at Z = 


(21b) 
ft, (21c) 


H = 


i 


n " 

■ — 

3 



Z(2h-Z) 

in 

12 



(4ft 2 + 2hZ - Z 2 



Z l 



3h-Z- — (20ft 3 - 5hZ 2 
20 v 



Z 2 

i-mlh 2 

12 



2 9 

-iVLh 2 
5 



(24a) 

(24b) 
Z 3 ) 
(24c) 
(24d) 



x 
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Equations (|24|) can also be obtained by expanding 
Eqs. ([H]) in powers of the small parameter a and keep- 
ing all terms of order a 4 . The terms proportional to fl 2 
originate from the evolution of U in t. These terms are 
small corrections, but they must be retained since they 
govern the solution of the averaged problem. 

In the fl 3> 1 case the amplitudes of the pulsations are 
[see Eqs. 



q = Q(h-Z), I = ihx, K = -ihxxZ, H = - (hh x ) x ■ 

(25) 

Equations (|25p constitute the solution to the vibration 
problem for an inviscid film. The longitudinal compo- 
nent of the velocity /, is uniform across the liquid layer; 
it is determined by the pulsations of the pressure gradi- 
ent only. The transversal component, K, is the linear 
function of Z; it vanishes at values of X corresponding 
to the extremum of I. Note that H is real, i.e the surface 
deformation is in phase or in antiphase with the vertical 
motion of the substrate. 

Of course, I given by Eq. (|25p is inconsistent with the 
no-slip condition (|21b|) . In order to vanish / at the rigid 
wall one has to take the boundary layer into account. 
Introducing the "fast" coordinate 77 = Z/yQ near the 
wall we arrive at 



0, u {i) - jw = _n- 

= at 77 = 0, 



K {i) 



0, 



(26) 
(27) 



ih x , q 



(i) 



q(Z = 0) = Qh at 77 -> 00. (28) 



The solution of this problem is well-known (see, for ex- 
ample, Ref. Hi): 



q(Z = 0) = Qh, 
ih x [1 - exp (-/?r?)] , 



0. 



(29a) 
(29b) 



where /? = Vi = (1 + i)/y/2. \I {i) {rf)\ increases from 
zero at the rigid wall to the maximal value 1.069|ftx| at 
77 as 2.284 and then decays to |/| = \hx\- Of course, the 
solution (J25JI matches the solution at O > 1, Z — > 0. 

Examples of the distribution of I r and Ii across the 
layer are given in Fig.[^a) and Fig.^b), respectively, for 
different values of il. (Here and below we use the sub- 
scripts "r" and "i" for the real and the imaginary parts 
of complex variables.) Since / is proportional to hx, the 
value of the latter derivative only rescales the longitudi- 
nal velocity. Thus we set hx = 1 for these sketches and 
for Fig. 1 

One can see immediately that at ft < il c w 19.74 the 
intensity of the oscillations is maximal at the free surface 
and monotonously decreases to the rigid wall. At higher 
values of Q there exists a maximum in the inner part of 
the layer. With increase of fl the Z-coordinate of this 
maximum tends to zero. At 17 = 100 the velocity profile 
agrees well with the asymptotic formula ([29)1 . 

The dependence of the maximal amplitude of the pul- 
sation velocity and of \I(Z — h)\ on the frequency of 
the vibration is shown in Fig. [3] by the solid and dashed 




(b) 



FIG. 2: Amplitude / of the longitudinal component of the 
oscillatory (pulsatile) velocity as function of Z/h. Real (a) 
and imaginary (b) parts. Lines 1-4 correspond to Q — 
0.5, 1, 10, 100. 




FIG. 3: Maximum longitudinal velocity of pulsations I m = 
max(j/|) as function of Q (solid line); the amplitude of this 
velocity at the surface (dashed line). 



lines, respectively. At large tt the solid line reaches the 
asymptotical value 1.069, the dashed line tends to unity. 

To make more clear the behavior of the pulsation veloc- 
ity we show in Figs. [4] and [5] the isolines of the pulsation 
streamfunction at four progressive time moments. The 
strcamfunction ^ is defined as 



* = ihx 



Z 



sina(/i — Z) — sin aft. 



a cos aft. 



(30) 



so that I — if?z, K = — $>x- Again for illustration pur- 
poses only we set ft = 1 + acoskX with a = 0.1, k = 1 
in these figures. 

We point out that these figures present the streamfunc- 
tions of the time-periodic motion, i.e. one must not be 
concerned that the isolines are open curves. Each fluid 
particle oscillates near its mean position with the small 
(on the scale of the figure) amplitude, and its instanta- 
neous pulsation velocity is tangential to the momentary 
isoline at the point. 

It must be emphasized that within the framework of 
the longwave approximation there is no need for the sta- 
bility analysis of the solution (f2"2"]l . Indeed, the leading 
order of Eqs. (Ti"4")) [given by Eqs. (TT5)) ] is the linear prob- 
lem. As was mentioned above, the homogeneous problem 



FIG. 4: Streamfunction \& of pulsations for fl — 1 at r = 
0, 7r/4, tt/2, 3tt/4 (top-to-bottom). 



FIG. 5: Streamfunction ^ of pulsations for = 10 at r = 
0,7r/4, 7r/2,37r/4 (top-to-bottom). 



originating from Eqs. (fl"9|) obviously has only the decay- 
ing solutions: at finite Q, the perturbations decay due 
to the viscosity. Thus, stability of the oscillatory mo- 
tion is evident except for the limitin g ca se fl 3> 1, which 
has been studied in detail [H, l34l. I36L However, 
stability with respect to perturbations of finite or small 
wavelength must be addressed. We briefly discuss this 
issue in Sec. IVIII Also, in the Appendix [2] we discuss 
the reduction of the flow ([22]) to the oscillatory Poiscuillc 
flow. The analysis in Sec. IVIII is partially based on this 
result. 



IV. AVERAGED MOTION 

In this section we solve the leading order of the aver- 
aged problem (fl6|) : 

Pz = ~G , W z = -U x , (31a) 

Uzz = Px + (UU X + WU Z ), (31b) 

U = W = at Z = 0, (31c) 

p = <p(h) - Chxx + bn 2 (hcosT), 

h T = W -Uhx + (hWz -Uh x ), 

U z = -(U zz h) at Z = h. (31d) 

Note that the averaged term in the boundary condition 
for the pressure results when pz in the corresponding 
boundary condition (|16d[) is replaced by — bVL 2 cos r [see 
Eqs. CEU)]. 

Substituting the forms (|20|) . accounting for the obvious 
equality 

(Re {Be") Re (De iT )) = ^Re (BD*) (32) 
for the calculation of averages in Eqs. (f3"Tj) . and noting 
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that 

(Uhx - Wzh) = (Uhx + U x h) = (Uh)x, 
we obtain (the overbars are omitted): 



Pz 

Uzz 
U 

P 

JlT 

U z 



-Go, W z 
Px 

W = at Z 
<P - Chxx H 



-U x , 
h 2 n 2 Re(I*I x 

0, 



K*I Z ) ■ 



(33) 

(34a) 
(34b) 
(34c) 



-^b 2 ilRe(I*H) x 



-Uh x + W 
- -b 2 URe {1 * ZZ H) at Z 



(34d) 



The evolutionary equation for h can be rewritten in the 
form 

f h 1 

h T = -d x / UdZ - -b 2 QRe (I*H) X at Z = h. (35) 
Jo 2 

Analytical integration of this set of equations is per- 
formed in Appendix [Bj It results in the following non- 
linear equation for h: 



~h 3 n x - Wq 



x 



1 



n = p(z = o) 

= <t>{h)-Chxx 
Q = Q 1 { 1 )h 2 h\ + 

where 7 = \/2£lh and 

2 sinh 7 sin 7 — 7 (sinh 7 cos 7 + sin 7 cosh 7) 
3 — ; — 72 ^* 



(36a) 



G h + -b z n z ReH. (36b) 
2 {j)h 3 h x hxx, (36c) 



Qi 



7 2 (cos 7 + cosh 7) 



■7a) 



1 1 1 (sinh 7 -sin 7) -37 (cosh 7- cos 7I 

— -- H 37 ; 7-^: b^Dj 

3 7° (cos 7 + cosn 7) 

Also note that the function f(ah) [see Eq. 
expressed in terms of 7: 

sinh 7 + sin 7 



f r (ah) 
fi(ah) 



1 



7 (cos 7 + cosh 7) 
sinh 7 — sin 7 



] can be 

(38a) 
(38b) 



7 (cos 7 + cosh 7) 

The dependence of functions f r ^(ah) and Q1.2 on the 
parameter 7 is given in Figs. [DJa) and^b), respectively. 
One can see that < f r (ah) < 1 and the coefficient Q2 
is negative for all values of 7. 

One can immediately see that along with the regular 
contributions due to the wetting potential, surface ten- 
sion and gravity, the expression for II contains the non- 
linear contribution due to the vibration. The nonlinear 
Q-term is entirely due to the vibration. 

The obtained equations of slow motion allow simpli- 
fication in the limits of low and high frequency f2. For 
these limiting cases, considered next, the solution to the 
boundary value problem (|34[) is not very cumbersome and 
can be presented in detail. 
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FIG. 6: Coefficients in the amplitude equation (|36[) as func- 
tions of 7 = y/2Qh: (a): Real and imaginary part of f(ah), 
(b): Qi,2(7). 



V. FREQUENCY-BASED ANALYSIS OF EQS. 
(HU) AND (|36l> 

A. Low vibration frequency, Q 1 

We look for the solution of the averaged problem (|34|) 
in the form of the power series in f2, using the corre- 
sponding solution (|24[) of the pulsatile problem. 

It is clear that W can be easily expressed via U by 
means of the continuity equation, the last relation in 
Eqs. (|34a|) . However, the corresponding expression in- 
volves different combinations of derivatives of h and is 
difficult to understand. Thus, for the sake of brevity, 
here we present only the longitudinal component of the 
averaged velocity: 



U = ^Z(Z-2h)U x 



(bn 2 



720 



-Z [75h 3 x hZ 3 



n 



-h x hxx (Z 5 - 6hZ 4 + 15h 2 Z 3 - 48/i 5 )](39a) 

(bn 2 ^ 2 



(j){h) - Ch 



xx 



G h 
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-{h b h x )x. (39b) 



(Also note that the pressure p = II — GqZ.) The terms 
proportional to f2 4 are retained in this solution, while the 
higher order terms have been omitted. The first term in 
Eq. (|39a[) is conventional, while the second one represents 
the impact of the vibration (as well as the last term in 
II). 

Substituting Eqs. ([39)) in Eq. (f35|) . we obtain the evo- 
lutionary equation for h: 

h T = i (h 3 n x ) x +(bn 2 ) 2 ^h 7 h x h X x - y^ 6 ^ 

(40) 

Equation (|4U|) can be obtained also from Eq. (|36p by 
noticing that at small 7 the coefficients Qi,2 (and f r ) 
are proportional to 7 4 : 



f r {ah) w Q\ 



-T 4 

30 ' ' 



63 ' - 



(41) 



Equation (|4H)l makes it clear that it is necessary to pro- 
vide large vibration amplitude in order to gain the fi- 
nite impact of the vibration on the dynamics of the film 
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height. In other words, the rescaled acceleration bfl 2 has 
to be finite. Moreover, we show in Sec. I VIII that such 
values of the amplitude do not cause the parametric in- 
stability. 

It is important to recognize that one needs only 7 <C 1 
in order to obtain Eq. PU)) . However, this limit can be 
reached not only for small f2, but also for small local 
values of h. This means that Eq. (T4"0|) can be used (in- 
dependently of the value of SI) near rupture, i.e. in the 
close vicinity of the point X r , such that h(X r ,T) — > 0. 
Note that the general Eq. (136)) should be applied far away 
from this point. 

However, since the vibration terms in Eq. (I40[) are pro- 
portional to h 9 , they are negligible at small h in compar- 
ison with the term originating from the surface tension 
(which is of order h 4 ) and with the dominant van der 
Waals term. In other words, only the competition of the 
surface tension and the van der Waals interaction governs 
the behavior of the film near rupture, and the vibration 
does not provide a noticeable impact. (Of course, in the 
very close vicinity of the rupture point only the van der 
Waals interaction contributes to the film dynamics - see 
Refs. PH].) 



B. High vibration frequency, f2 3> 1 

Here we use the 'inviscid" solution for the pulsations 
given by Eqs. (|2l)|) . We do not consider the impact of the 
boundary layer, i.e. the solution given by Eq. (I29j) . for 
the obvious reason. It is known from Rayleigh [421 ] and 
Schlichting 22], that the boundary layer can produce an 
independent averaged motion. However, the intensity of 
this flow is rather small in comparison with the volumet- 
ric sources under consideration. Indeed, estimating the 
longitudinal component of the averaged velocity gener- 
ated at the external border of the boundary layer, one 
obtains Us oc & 2 f2 (|/| 2 ) x , whereas the dominant contri- 
bution in the Eq. (|34ap is proportional to & 2 fi 2 (|/| 2 ) ^ . 
A similar situation exists in many problems of thermal 
vibrational convection [l7j ]. 

Now one can see that in the set (|34|l 

ReH = H = - (hh x ) x , (42a) 
Re (I*I X + K*I Z ) = h x h X x = \ (h x ) x ■ (42b) 

The corresponding terms are proportional to f2 2 . Other 
averaged terms in Eqs. (j34|) are proportional to £1, and 
thus they can be safely neglected. 
The calculation gives 



p = U-GZ, 

n = 4>(h)-Chxx+G h 



U = -Z(Z-2h) (n 



.202 



b 2 n 



(43a) 
(hhx)x, (43b) 



6 2 2 



h\ 



(43c) 



which in view of Eq. f|35[) results in the following equation 
governing the evolution of the thin film thickness: 



h T = ^{h 3 [^h)~Chxx + G Q h 



202 



b 2 n 



(2hh 



XX 



h 2 x ) 



X) X 



Since for 7> 1 

f r {ah) Qx = e.s.t., Q 2 



(44) 



(45) 



[36]) in the high 
denotes the ex- 



Eq. (|44[) also follows directly from Eq 
frequency approximation. Here "e.s.t.' 
ponentially small term. 

Obviously, at finite b and large the vibration deter- 
mines the film dynamics at large. One has to assume 
that b <C 1 to retain the competition of the surface ten- 
sion and gravity in the evolution of the film. Such small 
amplitude, high frequency vibration has been subject of 
many papers |l9l . I20I [32] ] . Equation (|4"4")) coincides with 
the equation derived, in the high frequency approxima- 
tion, by Lapuerta et al. [lj| (see also Ref. [2(|). For- 
mally, we obtain that the same equation remains valid at 
larger amplitudes, but in fact finite values of b cannot be 
reached for 3> 1 due to the instability of the pulsatile 
motion (see Ref. [37| and Se c.lVIIj) . 

Finally, we recall from Sec. Ill Bl that the high frequency 
limit is equivalent to the approximation of thick fluid 
layer. Thus the Eq. (|3S|). strictly speaking, can be applied 
only with <f> — 0. 



VI. 3D CASE 

In this section we generalize the theory to the 3D case. 
The starting point is the 3D analogue of Eqs. (jHJ) and ([5]): 



W z = 

nu T = 



e 2 (V 2 U 



U 7 



u 



(46a) 
(46b) 

Pz = -Go-Ml 2 cosT + e 2 (iy zz -fW T ) (46c) 
= 0, W = at Z = 0, (47a) 



-V-U, 
- Vp + Uzz 
-e 2 (U ■ VU + WUz 



Uz = -e 2 



\7W - 2V/iV • U - Vh ■ ( VU + VU T 



flh T 
P 



e 2 (W - h T - U • Vh) . 
6 - CV 2 h 



-C(\7hf\7 2 h + 2W Z 



at Z 



(47b) 



x 



Here V = (dx,dy,0) is the 2D gradient, U is the 
projection of the velocity onto the X — Y plane, i.e. 
v = eU + e 2 We z , (VU) is second-order tensor [i.e., 
(VU) j7 = dUi/dxj, j,l — 1,2], and other notations are 
unchanged. 
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It is easy to see that in the leading order (e°) the only 
difference between the systems (JS| and © (2D case) and 
(46]) and ([47]) (3D case) is the replacement of U by U 
and dx by V. Less evident is that only same changes 
are warranted in the solution as well. This can be easily 
checked by repeating the analysis quite similar to the one 
presented in Sees. iHlfVl Here we show the results only. 

The solution of the problem for the pulsations is [cf. 
Eqs. dH])]: 



Q(h-Z), I = i 



1 - 



cos a(h — Z) 



K 



H 



Z 



cos aft 
sina(/i — Z) — sin aft. 



V/i,(48a) 



a cos aft. 



V 2 ft 



1 



cosaZ a 



cos 2 aft 
-V • [ft/(aft)Vft] , 



(48b) 
(48c) 



while the averaged dynamics of the free surface is gov- 
erned by the following equation: 



ft T = V • 



ft 3 vn 



b 2 n 2 



Q(h) 



n = p(z = o) 

h 2 Q 2 

= 0(ft)-CV 2 ft + G o ft+- ReH. 



(49a) 



(49b) 



Q = Qi(7)ft 2 (Vft) 2 Vft + Q 2 (7)ft 3 V 2 ftVft, (49c) 

with Qi j2 and / given by Eqs. (j3"T]) - (58"]) , respectively. 
In the limiting cases Eq. (|4"9"]) simplifies as follows: 
(i) Low frequency, Q -C 1 



h T = -v • (ft 3 vn) 
+ (bn 2 ) 2 v 



2ft 7 



. 63 

n = 4>{h) - cv 2 ft + G ft- 

(ii) High frequency, f2 3> 1 



-V 2 ftVft 



ft 6 

— (Vft) 2 Vft 
15 V ' 



(,50) 



,2o4 



b 2 n 

15 



■V-(ft 5 Vft). (51) 



/i T = ly • <! ft a v 



n 



^<v.>* 



n = ^(ft) - cv 2 ft + G ft - 



b 2 fl 2 



(52) 

(ftVft) . (53) 



Omitting the van der Waals interaction potential </>(ft) 
one immediately reduces the last equation to Eqs. (3.11) 
and (3.12) in Ref. [H [see also Eq. (3.19) in Ref. at 
Ma = Q]. 



VII. STABILITY OF THE PULSATILE FLOW 
A. Analysis of the deformable mode 

In this section we analyze stability of the pulsatile mo- 
tion given by Eqs. (|22|) . It was mentioned already in 



Sec. IIIII that there is no any doubt about its stability 
within the framework of the longwave approximation, 
Eqs. (TTJ)) . However, the question whether this motion is 
stable with respect to perturbations with shorter wave- 
length is quite reasonable, especially in view of unusual 
scaling (jTTJ) : the amplitude of the vibration is large (of 
order e _1 ) and one can expect the emergence of the para- 
metric instability. 

Therefore, we return to the governing equations (3]) 
and (0| and study the stability of the "base state" : 



Po 

Wq 



= e -1 6fiRe [<?(ex, *) exp (ifi*)] + 0(1), 
= bQRe[I(ex,z)exp(int)]+0(e), 
= O{e),h Q = h (ex 1 e 2 t) + O(e). 



(54) 
(55) 
(56) 



[See Eqs. (|12bj) and (13]) for scalings and Eqs. (20]) for 
the pulsation velocity and pressure.] We restrict analysis 
to the 2D base state since the stability problem for the 
3D base state [fto = ho(ex, ey, e 2 t)} admits reduction to 
the one with the 2D base state. This will be shown be- 
low. Introducing small perturbations and linearizing the 
problem near the base state we obtain: 



u x + w z = 0, 



u t + u Q u x + wu 0z 

w t + U W X 



- Px + V 2 w, 
-p z + V 2 w, 



(57a) 
(57b) 
(57c) 

(58a) 



u = 0, w = at z = 0, 
& = w-unix, u z +w x = -u 0zz £, 
p = [e^Wcosfti + G Q + <p'(ho)] £ 

-e~ 2 C£ xx - 2w z = at z = ft , (58b) 

where £ is the perturbation to the surface deflection, and 
the obvious notations {p, u, w} are used for other pertur- 
bations. Again, we will show below that the 3D pertur- 
bations do not warrant the consideration. 

Note that for this analysis one can safely neglect vari- 
ation of fto on the time scale T — e 2 t and on the length 
scale X — ex. Thus the unperturbed surface is assumed 
to be locally flat. Same approximation is also appro- 
priate for the velocity components, thus we consider the 
plane-parallel flow with vanishing transversal component 
wo and the longitudinal component nearly constant in 
x. The obtained problem is quite similar to the prob- 
lem gov erning the Faraday instability (see, for example, 
Ref. j37|). The only difference is the presence of the base 
flow u in Eqs. (57]) and ([55]) . 

Presenting all fields in the form of normal pertur- 
bations (u,w,p,(,) = (u(t, z),w{t, z),p(t, z),i(t, z)J e lkx , 

where k is the wavenumber of the perturbation, we ob- 
tain the following problem for the amplitudes: 



iku + w z = 0, 
Ut + ikuou + wuq z = —ikp + u z 
Wt + iku w = -p z + w zl 



(59a) 
- k 2 u, (59b) 
k 2 w, (59c) 
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u = 0, w = at z = 0, (60a) 
| t = i(> - iku £, u z + ikw = -u 0zz i, 
p = [e- 1 bQ 2 coant + G + cf>'(ho)]i 

+e- 2 Ck 2 i - 2w z = at z = h . (60b) 

This problem contains terms of different orders with 
respect to small parameter e, which simplifies the analy- 
sis. Indeed, the term e~ 2 C^ xx prevails in the boundary 
condition and produces the stabilizing effect, since it sup- 
presses deviations of the surface. Instability may occur 
when this term is comparable to the potentially desta- 
bilizing term e~ 1 bQ 2 cosili in the same boundary condi- 
tion. This is only possible for long waves with k = yfiK. 
It can be shown easily that for other fc the perturbations 
decay, except for the case considered in Sec. IVH B[ 

Choosing the following scalings for the perturbations 
(see Ref. [371]): 



p = P, u = \/~eJJ , w = eW, | = eS, 
we obtain in the leading order 



(61) 



-iKP + U zz , iKU + W z = 0, P z = 0, (62a) 
0, w = at z = 0, (62b) 
W, U z = 0, 

bn 2 cos SltE + CK 2 E = at z = h . (62c) 



U t = 

u = 

P = 



It can be seen that in Eqs. (|59|) and (1601) the terms con- 
taining base flow uq are of low order, and thus they 
dropped out of Eqs. I]62p . This means that within the 
framework of scaling (|6ip Faraday instability prevails 
over the instability due to shear flow (see Sec. I VII Bl for 
the opposite case). 

Next, it is evident that there is no preferential direc- 
tion in the x — y plane for the problem (|62p . This allows 
the reduction of the stability with respect to 3D pertur- 
bations to 2D problem under consideration: one needs 
only to choose the x-axis in the direction of perturba- 
tions' wavevector. Moreover, since the velocity of the 
base state also dropped out from the leading order of the 
stability problem, the symmetry properties of the base 
state are inessential. Therefore, even for 3D base state 
the stability is governed by the same Eqs. ((62)) . 

Thus, we showed that the problem under consideration 
completely reduces to the analysis of Faraday instability. 
Such analysis was performed in detail by Mancebo and 
Vega [37]. Equations (62) can be rewritten in the form 
of Eqs. (2.11)-(2.13) in their paper. Indeed, the case 
considered here corresponds to the case B.1.2 (Long-wave 
limit) in their paper. Introducing the parameters as in 
Ref. 13711 . we obtain: 



flh , clmv 



Ch 



3/2 



Go ,5/2 Vco -. 
1MV = £ -^ h < X 63 ) 



c 



Here the subscript "MV" is used to mark the parame- 
ters used by Mancebo and Vega [37j • Note that these 



parameters must be calculated using the local thickness 
H = HohQ instead of the mean value Hq, which causes 
the appearance of ho in Eqs. (|6"3"]) . 

The critical value of the acceleration aMv. A v as a 
function of ljmv is presented in Fig. 4 of Ref. (37|. Thus, 
the critical value of the amplitude is 



(64) 



C^mv^IK) 
Q 2 h 5 J 2 



where $(/i) is the function given in Ref. [33] (Fig. 4 
there). Due to Eq. (|63|) we are interested only in the 
line corresponding to Jmv = 0. 

At \i — > 0, $mv(aO ~ 8.5, then it decreases to approx- 
imately 5.5 at p « 6 and after that grows. At large ujmv 
the asymptotic formula 



holds. It means that 

b c (Qh f 2 C-^ 2 



1 at n 



(65) 



(66) 



Using Eqs. ([55]). (|4"0"]1. or (gU) one has to ensure that 
b < b c . Value of b c can be extracted from the men- 
tioned figure at given fl. Note that the stability condi- 
tion should be valid at any X, therefore the value of hg, 

which minimizes $ mv (^o)^'O^ 2 eacn time moment) 
must be substituted in Eq. (|64|) to determine whether 
the layer is stable or not. For most of cases, except for 
3 < Q < 8, this means that the maximum value of ho 
should be used. In the opposite case the minimum value 
of $mv(/j) (~ 5.5) and again the maximum of ho can be 
used in order to estimate b c . 

For the case of low frequency one can take the acceler- 
ation bfl 2 up to 



b c VL z « 8.5VCK 



-5/2 



(67) 



in Eqs. ([101) . I n the opposite limiting case, S! 3> 1, the 
parameter M7, which enters Eq. ([4"4"[) , should be small: 



M> b,\l = \f ^hl < 1. 



(68) 



The inequality ([6"8"[) holds in the limiting case 1 -C H <C 
e" 1 , where the Faraday instability is caused by the long- 
wave perturbations and the dissipation is negligible be- 
yond the boundary layers near the rigid wall and free 
surface. Similar to Eq. (l6"5[) the stability bound b c fl is 
small even for Q = (9(e _1 ); the perturbations with mod- 
erate wavelength are critical in this case. Thus, in this 
case 60 must also be small to prevent the parametric 
instability. 

For very large frequency, i.e. f2 3> e _1 , the volume 
dissipation prevails, which leads to another limitation on 
the amplitude. Indeed, it follows from Ref. [37j] that the 
inequality 



l'd-< A MV ( - J -^ r - 



(69) 
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(in dimensional form) or 



(70) 



(in dimensionless form) must hold true to prevent the 
Faraday instability. Here Amv(u) is the function ob- 
tained by Mancebo and Vega [371] , which has the follow- 
ing asymptotics: 

A M v (m) « (256^) 1/6 at /x < 1. (71) 
This gives the following critical velocity: 

/256fte 



\ C 2 



\ C 



(72) 



Multiplying this relation by we obtain: 

'3 

bcM = [ ^7- ) , (73) 



C 



i.e. the vibration is non-negligible in Eq. (|44"]) at high 
frequency only if f2 2 e 5 is finite or large. This provides, 
of course, the usual limitation on the amplitude and the 
frequency of the vibration in the case of inviscid pulsatile 
motion (see also Refs. 0, Ql, S3 ) . 

It must be emphasized that advective terms, being pro- 
portional to uq, are not important for short waves, since 
the advective term ikugu remains small in comparison 
with ut because of the dispersion relation u> 2 = Cak 3 at 
large k. 



B. Analysis of the nondeformable mode 

The above stability analysis deals with the capillary 
waves, which are based on the surface deviations. lib < 
b c , where b c is given by Eq. ((64)) , the perturbations of 
this type decay. 

However, there is also a mode, which corresponds to 
the nondeformable surface. Indeed, it is obvious that the 
finite-amplitude plane-parallel flow with the profile uq 
becomes unstable at certain intensity of the motion. Let 
us analyze this mode in detail. Setting £ = in Eqs. ([5"5]) 
and (]6TJ]) gives: 



ft, 



iku + w z 
ikuQii + wuqz 
w t + ikuQW 
u 

W 



0, 

—ikp 
~Pz 4 
0, w -- 
0, u z 



- V 2 w, 
V 2 w, 

= at z = 0, 
= at z = ho. 



(74a) 
(74b) 
(74c) 
(74d) 
(74e) 



This problem is one of stability for the periodic in time, 
plane-parallel flow. In view of Eqs. (|22|) and ([54]) this 
problem is characterized only by the parameters b, 
and "local" values of ho and hox- Due to well-known 
Squire theorem [43| there is no need in analysis of 3D 



perturbations - 2D one are critical. Moreover, the base 
state should not necessarily be 2D in the entire layer. It 
is sufficient that the flow is locally 2D at any point in the 
X — Y plane. 

Note that according to the analysis in Appendix [5] 
the problem ([71]) for the base flow ([5"4"]) . but posed on 
the interval < z < 2ho and with the no-slip con- 
dition at z = 2ho instead of Eq. (|74e|) . is the conven- 
tional stability problem of the oscillatory Poiseuille flow 
(i.e., the flow which arises due to the periodically os- 
cillating longitudinal pressure gradient P x = Pi cos fit) . 
The more general problem, where the pressure gradient 
equals Pn + Pi cos fii, has been investigated already, see 
Refs. [io, EH, H3] and references therein. Due to sym- 
metry the latter problem can be split into problems for 
"even" and for "odd" perturbations, meaning that w is 
even or odd function of the coordinate z — ho- For the 
odd mode both w and u z vanish at z = ho, which co- 
incides with the boundary condition (|74ep . Therefore, 
the problem under consideration is the particular case of 
the stability problem for the oscillatory Poiseuille flow. 
However, to the best of our knowledge there is no de- 
tailed analysis of the flow stability for the particular case 
we need, i.e. Po = 0. 

It is easy to consider two limiting cases, i.e. the high 
frequency and the low frequency. The first case, in view 
of Eq. ([23]) and ([2"5]) is reduced to the stability of a Stokes 
layer. The latter is known to be stable 0, |4f|. 

At il -> we have u oc bVL 2 from Eqs. (gD]) and ([2"4"]l . 
On the other hand in this case we can "freeze" the evo- 
lution of the flow and assume that 



bVt 2 



u 



h x z(2h - z). 



(75) 



This means that we assume the frequency so low that 
the perturbations either grow or decay before the flow in 
fixed point of the layer will change itself. 

In view of the above-mentioned symmetry properties, 
the problem ([74]) with u given by Eq. ([75]) is identical to 
the stability problem for the stationary Poiseuille flow in 
the entire layer (ho being the half of the layer thickness), 
but only for the odd perturbations. If we introduce the 
Reynolds number based on the velocity of the flow at 
z = ho, we obtain: 



Re 



bn 2 h x hl 



(76) 



[It follows from Eq. ([74]) that the viscosity is equal to 
unity] 

It is known [46| that the critical value of the Reynolds 
number is 5772, i.e. the flow remains stable for 



Re < Re r = 5772. 



(77) 



Thus the flow forced by low frequency vibration is stable 
if 



bfl 2 < b cl n 2 = 



11542 
hxhl ' 



(78) 
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Of course, this limitation is less severe than Eq. (|67|) for 
any reasonable value of hx- The product hxhg should 
be maximized over the longitudinal coordinate X at each 
moment of time. Thus, the pulsatile flow is stable at low 
frequency up to finite value of the acceleration, i.e. the 
vibration can produce the finite impact in this limiting 
case. The complete analysis of the problem (l?4")) will 
be performed elsewhere. This research will provide the 
threshold value of the amplitude 



Re c (Qhl) = b cl Q\I{ah )\h h 



x- 



(79) 



(We again define the Reynolds number via the velocity 
at the center of the layer, z = ho.) Recall that due to the 
minimization of b c i, the maximal value of \I{otho)\hohx 
should be used in Eq. ([75]). 

To conclude, we showed in this section that there exists 
an upper bound b = b c below which the pulsatile flow 
is stable. The further analysis as well as the results of 
Sees. IIIlllVIl are based on the assumption b < b c . 



VIII. 



EVOLUTION OF THE PERTURBATIONS 
TO THE FLAT LAYER 

A. Linear stability analysis 



The amplitude equation (|3l)|) has the obvious solution 
ho = 1 corresponding to the equilibrium state. It follows 
from Eqs. ^ and jM]) that 

ho = 1, po = G (l -Z) + 0(1), q = 0(1 - Zj80a) 
U Q = Wo = Io = K = H = 0. (80b) 

Thus in the reference frame of the substrate the fluid is 
motionless - vibration only adds the oscillatory compo- 
nent to the pressure field. 

Let in Eq. (j3"6")l h = l+£, where £ is small perturbation. 
Linearizing with respect to £ we obtain: 



= ~ <!;<<'( :i ) + G,i- 



c 



b 2 n 2 f r (a) 



XX 



XX 

(81) 

where a 2 — —iCl and f(y) is given by Eq. (j23|) . 

The linear stability of the layer without accounting for 
van der Waals attraction, and for the opposite direction 
of the gravity field was studied by Lapuerta et al. [l9[ . In 
this setup the Rayleigh- Taylor instability emerges. The 
equations governing the dynamics of small perturbations 
derived in Ref. [lj| [see Eqs. (2.35) there] can be obtained 
from Eq. (fgTj) . 

However, even in this case there is an important dif- 
ference in the interpretation of the results. Lapuerta et 
al. [l|| address the case of finite vibration amplitude B 
and consequently, they find that the influence of the vi- 
bration is small. To gain a finite impact of the vibration 
they proceed to the detailed analysis in the limit of high 
vibration frequencies, f2 3> 1. Conversely, in this paper 
large vibration amplitude is considered and we show that 



the Eq. ([81]) remains valid even in this case. Thus, we 
extend the domain of applicability of the results obtained 
by Lapuerta et al. by showing that the finite impact of 
the vibration is possible even at moderate frequencies, 
which is most important for thin films. 

The typical stability curves are shown in Ref. [l9( : Fig- 
ure 4 there presents the dependence of the dimensionless 
amplitude of the vibration Al on the dimensionless fre- 
quency u>l for different values of the parameter a l , which 
is proportional to the surface tension. The results of our 
linear stability analysis (for cf> — 0) can be extracted from 
this figure when the following substitutions are made: 



A, 



bn 2 



3A-G 



, LU L — £1, OIL — 



k 2 C 



3A - Go 



(82) 



[Recall that 4>(h) = Ah~ 3 .] However, to avoid the recal- 
culation we present the results of stability analysis below. 

Seeking the solution in the form of a plane wave £ — 
£o exp (—AT + ikX), where A is the decay rate and k is 
the wavenumber, we obtain 



A = \k 2 U'(l) + Go + k 2 



C+h 2 n 2 f r (a) 



The stability criteria (A = 0) is 

1 



0'(1) + G + 



C+-b 2 Q 2 .f r (a) 



k 2 =0. 



(83) 



(84) 



Thus the vibration and the surface tension do not damp 
the longwave instability: the perturbations with small k 
grow at cf>' (1) + Gq < 0. However, in confined cavities 
with large aspect ratios the spectrum of the wavenum- 
bers is discrete and bounded from below. Therefore the 
impact of the vibration and the surface tension becomes 
determinative in this case. 

From Eq. (|84p one can see that the critical value of the 
wavenumber, k c (i.e the value that corresponds to van- 
ishing growth rate of the perturbation) becomes smaller 
due to vibration: 



#(1) + G 
C +lb 2 n 2 f r (a)' 



(85) 



Again we note that f r grows monotonically from zero 
(at Q — ► 0) to unity (at £1 — ► oo), see Fig. EJa). Thus 
the vibration leads to the stabilization of the thin film. 
This stabilization effect is obviously augmented with the 
increase of the frequency, even when this increase is ac- 
companied by the decrease of b to keep fixed the power 
of the vibration, b 2 fl 2 . 

At large Eq. |HJ reduces to 



(/>'(!) + Go + [C + 



b 2 ^ 2 



k = 0. 



(86) 



A similar equation (without the first term and with neg- 
ative Go) is the well-known result on the suppression of 
the Rayleigh- Taylor instability 0, [H, 13] • 



14 




FIG. 7: Linear decay rate of the perturbation, A [Eq. (|87|l ] as function of the wavenumber k. Go = 3.33 ■ 10~ 4 . (a): 0. — 
0.2, V = 0, 50, 200 (dashed, solid, and dash-dotted lines, respectively); (b): Q = 2; (c): Q — 20. For (b) and (c) dashed, solid, 
and dash-dotted lines correspond to V = 0, 1, 5, respectively. 




FIG. 8: Critical value of the wavenumber k c as function of 
Vo = V{1~ 2 and the frequency Q. Go = 3.33 ■ 10" 4 . 



Introduction of the rescaled wavenumber k = k\JC '/3A 
and the growth rate A = AC/ (3A 2 ) results in the follow- 
ing expressions for A: 



A = fc 2 |g — 1 + k 2 [1 
and for the critical wavenumber: 

Here 



Vf r (a)]}, 



(87) 



1-Go 
l + Vf r {a)' 



r Go y _bW 



b 2 q 2 

2Ca 



(89) 



are the rescaled Galileo number and the vibrational pa- 
rameter, respectively. The Galileo number usually is 



quite small for thin films. We take Go = 3.3 • 10 4 . 
This corresponds to A' = 6ir ■ 10 -21 J (value for water) 
and Hq = 1000A Dependence of \(k) is shown in Fig. [7] 
for various values of and V. Clearly, k c decreases with 
V. Moreover, as it is shown in Fig. [51 k c also decreases 
with growth of Q,. Thus the vibration stabilizes the film. 
And, stabilization is more pronounced for larger values 

of a 



B. Weakly nonlinear analysis 

Let consider the behavior of perturbations near the 
stability threshold determined in the previous subsection. 
For this purpose it is convenient to rescale the time and 
the coordinate as follows: 



(90) 



Substituting these relations into Eqs. (|36|) one can obtain 



hf 



n 

Q 



X 



- h^x + G h + VReH, 
Q 1 { 1 )h 2 h%+Q 2 { 1 )h z h jl h^, 



(91a) 

(91b) 
(91c) 



where Qi,2 are given by Eqs. (|37j) and the rescaled pa- 
rameters defined by Eq. (|59")) are used. Below we use only 
the rescaled coordinate and the time. Thus we omit the 
tildes above X and T. 

In Sec. lVIlOl it has been shown that the growth rate is 
real at the stability threshold. Consequently, the branch- 
ing solution is stationary, i.e. one can omit the left- 
hand side of the amplitude equation to study direction of 
branching only. This also allows us to integrate Eq. (f9"Tj) 
once. Besides, seeking the solution with fixed wavenum- 
ber k we introduce the variable ( = kX. The surface 
deflection h now is a 2-7r-periodic function of C 5 an d it 
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solves the following equation: 



1 

3P 



k 2 h C c + Goh + Vk 2 ReH) -Wk 2 Q = const. 



(92) 

Note that X must be replaced with £ in the expressions 
for Qi i2 and H. 

Next, we expand h near the base solution ([80|) as fol- 
lows: 



which obviously coincides with Eq. (|88|) . Recall that we 
are interested in the case Gq < 1 which holds true for the 
thin film. 

The equation at the second order in 5 has the form: 



Lh 2 = 2(h{)' + VkfiFi (hitif + 3h[h'l] 
-WklQzKh'l. 



(98) 



h = 1 + 5hi + S 2 h 2 



(93) 



The wavenumber is assumed close to the critical value 
fc , i.e. 



fc = fco + <5 2 &2 



(94) 



Substituting these expansions into Eq. (|92[) , we arrive at 
order zero to 

Lhi = (1 - GoK + k 2 Q [1 + Vf r (a)} h'{' = 0, (95) 

where the prime denotes the derivative with respect to 
£. Solution of this equation is 

hi = a cos (. (96) 

The wavenumber is given by the expression 



P - 1 ~ G ° 



l + Vf r (a)' 



(97) 



The solution of this equation is 



h 2 = Cia 2 cos2C, 



Hereafter 



1 



3(1 -Go) 



h=l 



dhl 



(99) 



^= , F = hf r (ah). 



h=l 



At the third order we need only the solvability condi- 
tion. (We do not present here the corresponding equa- 
tion for /13.) This condition couples the correction to the 
wavenumber k 2 to the amplitude of the perturbation a, 
as follows: 



k a 2 



2 1 - G, 



- - 2d + Vk 2 I -Qi + 3Q 2 G! - -Q2 ~ 2 ClFl ~ g F2 



(101) 



Our numerical simulations show that the expression in 
the square brackets is always positive. Thus fc 2 > 0, i.e. 
a small amplitude solution exists at k > fco- This means 
that the subcritical bifurcation takes place. 

To summarize, the solution emerging at k = ko is un- 
stable and a finite amplitude excitation (probably leading 
to rupture) is expected. 

C. Nonlinear evolution of perturbations 

In this Section we analyze the finite-amplitude deflec- 
tions of the free surface, starting from stationary solu- 
tions h s of the amplitude equation. For this purpose 
we look for the periodic (with respect to £) solutions of 
Eq. (HU). 

Due to evident symmetry properties of the solution 
one can integrate Eq. over half of the period, van- 
ishing all the odd derivatives at £ = 0,7r. This means 



that const. = in Eq. (|9"2")l . We also have two boundary 
conditions: 

h' 8 = at C = 0, 7T. (102) 

Besides, the stationary solution conserves the liquid vol- 
ume: 




This provides the third boundary condition for the third 
order ODE, Eq. (f9"2"]) . completing the problem statement. 

The shooting method is applied to numerically inte- 
grate this boundary value problem. Some results are 
presented in Figs. [H] and [TO] It is clearly seen that 
only the subcritical (and, consequently, unstable) solu- 
tion branch is present, i.e. there is no bifurcation except 
the inverse pitchfork bifurcation studied in Sec. IVIII Bl 
Therefore, the lower branches in Fig. [TO] are the bound- 
aries of domains of attraction: the initial perturbation 
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FIG. 10: Minimal (lower branch) and maximal (upper branch) thickness of the film at the corresponding unstable station- 
ary solution. Go = 3.33 • 10~ 4 . (a): Q — 0.02, dashed, solid, and dash-dotted lines correspond to V = 0, 800, 5000 
(fco = 0.99998, 0.9793, 0.8885); (b): Q, = 0.2, dashed, solid, and dash-dotted lines correspond to V = 0, 50, 200 
(k = 0.99998, 0.8891, 0.6968); (c): fi = 2, dashed, solid, and dash-dotted lines correspond to V = 0, 8, 50 (fe = 
0.99998, 0.5287, 0.2418). Arrows show boundaries of domains of attraction. 



1 -i 

Kun 




5 t 10 

FIG. 11: Evolution of the minimum film thickness, as given 
by direct numerical computation of Eq. (|36[) . Unstable sta- 
tionary solutions with small perturbation of the free surface 
Sh were chosen for the initial condition. Go = 3.33- 10~ 4 , k = 
1.1, n = 0.02, V = 800. 5 h = 0.009, 0.002, 0, -0.002, -0.01 
(dash-double-dotted, dash-dotted, solid, dotted, and dashed 
fines, respectively). 

with min[/i((,T = 0)] > h m decays and the free surface 
becomes flat, while in the opposite case the free surface 
is attracted to the solid, which leads to rupture. These 
predictions are in good agreement with the results of the 
numerical simulation of Eqs. ([§1"]) , as shown in Fig. QT] 
The upper branches in Figs. [TO] are of importance for 



the stability of the oscillatory flow (see Sec. IVII[) . In- 
deed, Eq. ([6~4"]) requires the maximal value of the surface 
deviation ho, which can be extracted from Fig. 1101 

It is important to point out that the increase of the 
vibration amplitude amplifies the film stability: at fixed 
k and with Q increasing, the initial deviation of the flat 
surface decays at larger values of V. The only exception 
is provided by il = 2. In this case h m is nonmonotonic 
function of V starting from certain value of k. Also it 
is interesting that the maximal deviation of the surface 
decreases with growth of V: the surface tends to become 
flat as it is shown in Fig. [Sj In some sense this tend ency 
is reminiscent of the known averaged behavior [18|, [47| : 
the free surface/interface tries to orient normally to the 
vibration axis. But it is necessary to keep in mind that 
the surfaces shown in Fig. [5] correspond to unstable states. 



IX. SUMMARY 

In this paper the impacts of the vertical vibration on 
the dynamics of the thin liquid film are analyzed. The set 
of equations governing the averaged dynamics of the fluid 
flow and the nonlinear, fourth-order amplitude equation 
([56]) [or (|49p in 3D case] describing the averaged evolu- 
tion of the film thickness are obtained in the lubrication 
approximation. 
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We use the paradoxical scaling (at least at the first 
glance), assuming (i) the vibration period is compara- 
ble to the characteristic time of the momentum relax- 
ation across the layer and (ii) the vibration amplitude is 
large in comparison with the mean layer thickness. The 
first condition (termed "moderate frequency" or "finite 
frequency") allows us to consider ultrathin liquid layers 
within the framework of the averaging method. The sec- 
ond condition warrants that the impact of the vibration 
is not vanishingly small. Using the results from Refs. 
[37I [iol l4lj we prove that these assumptions do not nec- 
essarily lead to the parametric instability. 

Indeed, we show that the stability problem for the pul- 
satile flow can be separated into the problems for the de- 
formable and non-deformable modes. The former prob- 
lem reduces to the analysis of Faraday instability, while 
the latter problem reduces to the stability analysis of the 
oscillatory Poiseuille flow. Analyzing well known results 
obtained for these two problems we deduce that there 
exists a window of stability: below certain threshold in- 
tensity of the vibration [Eqs. and ([75)) ] the pulsatile 
flow is stable, while the averaged effects are well pro- 
nounced. 

The analyses of the averaged dynamics of the thin film 
demonstrate the strong stabilizing impact of the vibra- 
tion. First, the vibration damps the short- wavelength 
instability. In other words, it decreases the critical cut- 
off wavenumber k c , such that instability occurs at k < k c 
only [see Eqs. ([Ml), © and Figs. [JJ and [S]. In this sense 
the vibration acts in a way similar to the surface tension. 
Therefore, in order to prevent a longwave instability one 
can use a cavity of horizontal size L < L c ~ k" 1 , which 
is larger in presence of the vibration. Second, the vi- 
bration augments the domain of attraction of the flat 
undeformed surface, i.e. larger initial surface deflections 
decay (see Fig. [10]) or, in other words, larger initial distor- 
tions of a flat surface are admissible without occurrence 
of dewetting. 

Thus, the vertical vibration of moderate frequency 
is the effective method of control of the thin film in- 
stability. This is especially important since the stan- 
dard (high-frequency) approximation cannot be applied 
to thin films. 
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APPENDIX A: REDUCTION OF PULSATILE 
FLOW TO THE OSCILLATORY POISEUILLE 
FLOW 

In this section we show that the pulsatile flow given 
by Eqs. (|20[) and ([2"2"|) can be reduced to the well- known 
oscillatory Poiseuille flow. Such a transformation is use- 
ful for the stability analysis carried out in Sec. IVII Bl 
Throughout this section we omit the tildes over the ve- 
locity and the pressure field, with the understanding that 
only the oscillatory components are involved. 

In order to obtain the oscillatory Poiseuille flow we 
return to the unsealed coordinates x and z, setting locally 



h = ho + 0(e), hx = const + 0(e). 



(Al) 



These equations mean that the layer thickness changes 
slowly with x. Thus one can assume the constant thick- 
ness ho. In view of the scaling (fT2b|) and Eqs. {[3DJ 
and ([22]) the longitudinal pressure gradient is bSlhx cos r. 
Thus the pressure gradient is spatially uniform but oscil- 
lates in time. This means that some kind of the so-called 
oscillatory Poiseuille flow is under consideration. 

Neglecting hxx according to Eq. (|A1|) we arrive at the 
following expressions for the amplitudes of pulsations 



I = ih x 1 - 



cos aQiQ — z) 



cos aho 



(A2a) 



K = ih 



X' 



1 — cos aZ 
cos 2 aho 



, H = h\ tan 2 ah , (A2b) 



which along with the scalings (|12b[) and (IT3|) gives: 



cos a(ho — z) ' 
l[l : I e" 



u = bQhxRe - 

cos aho 
w = 0(e), h-h = 0(e). 



(A3) 
(A4) 



Thus this case corresponds to oscillatory ID flow in a 
locally flat layer under the spatially uniform longitudinal 
gradient of pressure. Note that the second relation in 
Eqs. (|A4|) justifies the first assumption in Eqs. (|A1[) . 

Separating the real and imaginary parts in Eq. (|A3p we 
arrive at the following expressions for the x-component 
of the pulsation velocity: 



-bQh 



x 



as{ho)a s (z) + b s (h )b s (z) 



b%(h ) 



a s (z) 



a s (ho) -r u s \ 
cos a r z cosh a r z, bs(z) = sina r z sinh cv r .z, 



as(h )b s (z) - a s (z)b s (h ) 



a 2 s (h Q ) 



b%(h ) 



COST 



(A5) 



Of 
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Such flow, but in a gap between two rigid boundaries, 
is well studied - see, for instance, the recent papers by 
Singer et al. [!(| and Straatman et al. [41], who ad- 
dress the stability problem for the flow. (Usually the 
modulated Poiseuille flow is considered, i.e. the pressure 
gradient oscillates about a non-zero mean value.) 

Due to symmetry the oscillatory Poiseuille flows in a 
layer with the free nondeformable surface at z = ho and in 
a layer with the upper rigid wall at z = 2ho are identical. 
Indeed, in the latter case in the plane of symmetry z = h a 
the "no-stress" condition u z — is obviously held. 

Thus, introducing variable ys = z — ho in Eq. (|A5|) 
we obtain the velocity profile of the oscillatory Poiseuille 
flow [cf. oscillatory part of Eq. (2.2) in Ref. |4(i[]. 



APPENDIX B: SOLUTION OF THE SET OF 
AVERAGED EQUATIONS 

Solution of the problem for the averaged fields can be 
represented as the sum of two solutions. The first one is 
conventional (see, for example, Ref. (HI); we mark it by 
the subscript "c": 



Pc 
U c 



= n-G z 

= <j)~Chxx+G Q (h-Z) + 

= \z{z-2h)n x , 



b 2 n 2 



1 



[z 2 (z-3h)n x ] 



x ' 



■ReH, (Bla) 
(Bib) 
(Blc) 



Note that this solution leads to the first term in the evo- 
lution equation for the film thickness, Eq. (|36[). This part 
of the solution coincides with Eqs. (2.47) in Ref. [ll| up 
to the averaged correction to the pressure field. 

The second part is the solution of the nonhomogeneous 
boundary value problem with the remaining vibration- 
generated terms at the right-hand sides. Using subscript 
"u" for this part of solution, we rewrite it in the following 
form: 

p v = 0, u v = h 2 n 2 u^ v \ w v = h 2 n 2 w^ v \ (B2) 



where and W^"' solve the following boundary value 
problem: 



Vf z = Re (PI X + K*I Z ) , W { z v) = -U%\ (B3a) 



Jjiv) = w (v) = q at z = Q) 

r(«) 



(B3b) 

U y z > = -Q^Re (i| z il) at Z = h. (B3c) 
This set of equations can be easily integrated. First, 

U ( z v) = -fi^Re [I* zz (h)H] - Re (I*I X + K*I Z ) dZ. 

J z 

(B4) 

Accounting for Eq. (f2"Tj) of the pulsatile motion and in- 
tegrating by parts, one can rewrite the last expression 
as 

U ( z ' = -h x ReH + Re (I*K) - d x \I\ 2 dZ. (B5) 

Jz 

After one more integration we arrive at the following so- 
lution: 

r z 

t/W = -Zh x ReH+ / Re(PK)dZ 
Jo 

-dx £ : d( J \I(X,0\ 2 d£. (B6) 



Evaluation of these integrals leads to the cumbersome 
formulas, which we do not present here. 

We also do not present the expression for W^ v \ as it is 
not needed in order to obtain the evolution equation for 
film thickness h. Indeed, this part of the solution results 
in the term 

r h 

-d x U {v) dZ -n^ReiPH} (B7) 
Jo 

at the right-hand side of such an equation. This term 
translates to the term Q(h) in Eq. ([55)) . 
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